clear
version 14
set more off
cap log close

mata: mata clear
mat drop _all

gl OUT 	"${LOCDRIVE}\Results"
gl TEMP	"${LOCDRIVE}\Tempdata"
gl FIG "${LOCDRIVE}\Figures"

loc r1_hh hhid
loc r1_weight hh_weight
loc r1_indiv sbmemno
loc r2_hh y2_hhid
loc r2_weight y2_weight
loc r2_indiv indidy2
loc r3_hh y3_hhid
loc r3_weight y3_weight
loc r3_indiv indidy3

loc ea ea
loc strata strataid
loc admin0 macroregion
loc admin1 region 		
loc admin2 district  	
loc admin3 ward

loc numeraire other

loc cre 1 /*include correlated random effects*/
loc crelev hh_c /*level of correlated random effects*/ 

loc servicelab "Services"
loc food1lab "SF"
loc food2lab "PNFFV"
loc food3lab "ASF"
loc food4lab "OF"
loc foodlab "Food"
loc goodlab "Goods"
loc otherlab "Other"

loc g1 food
loc g2 good
loc g3 service 

loc dobootstrap 1  /*Set to 1 to calculate the elasticity each iteration of the bootstrap, 0 otherwise if the dataset was already made previously*/

forvalues r = 1/3 {
	loc xc_dfl_r`r' = ${XC_DFL_R`r'}
	loc dfl_r`r' = ${DFL_R`r'}
	}
	
log using "${TEMP}\Analaysis_postest", replace


/*Set the main specification to be analyzed here*/
loc mainspec SL_food1good1service1_np5 


********************************************
** Define a program to recover parameters **
********************************************

cap program drop param_recover_easi

program define param_recover_easi, rclass
args pvector neq

noisily {
mat `pvector' = `pvector'
loc K = colsof(`pvector')
loc S = `neq' + 1 /*last (recovered) equation*/
loc Sminus1 = `neq' 

forvalues s = 1/`S' {
    cap mat drop `pvector'_s`s'
	}

forvalues s = 1/`Sminus1' {
    mat `pvector'_s`s' = `pvector'[1,1+(`K'/`Sminus1')*(`s'-1)..(`K'/`Sminus1')*`s']
*	mat list `pvector'_s`s'
	}

	
* Make a blank matrix for recovering parameters of the dropped equation 
loc colnames: colnames `pvector'_s`Sminus1'
mat `pvector'_s`S' = J(1,`K'/`Sminus1',.)
mat coleq `pvector'_s`S' = s`S'
mat colnames `pvector'_s`S' = `colnames'


*** Recover lost parameters

* first create colums for parameters to be recovered within equations
forvalues s = 1/`S' {
	foreach stub in np`Sminus1' ynp`Sminus1' np`Sminus1'z1 {
		loc index = colnumb(`pvector'_s`s',"`stub'")
		mat left = `pvector'_s`s'[1,1..`index']
		mat right = `pvector'_s`s'[1,`index'+1...]

		mat newparam = J(1,1,.) 
		mat coleq newparam = s`s'
		loc newparam_name = subinstr("`stub'","`Sminus1'","`S'",1)
		mat colnames newparam = `newparam_name'

		mat `pvector'_s`s' = left , newparam
		mat `pvector'_s`s' = `pvector'_s`s' , right 
		}
*	mat list `pvector'_s`s' 
	}


** Cross-equation sum to zero

loc sumlist = "$ylist $zlist $yzlist zs0 $zslist $creplist $crepylist"
display "`sumlist'"

foreach var in `sumlist' {
	loc index = colnumb(`pvector'_s1,"`var'")
	loc param_s`S' = 0 
	forvalues s = 1/`Sminus1' {
		loc param_s`S' = `param_s`S'' + `pvector'_s`s'[1,`index']
		}
	mat `pvector'_s`S'[1,`index'] = (0 - `param_s`S'')
	}
*mat list `pvector'_s`S'


** Summing to zero within s1 and s2 

foreach stub in np`S' ynp`S' np`S'z1 {
	loc param = 0
    forvalues s = 1/`Sminus1' {
	    loc index2 = colnumb(`pvector'_s`s',"`stub'")
		forvalues r = 1/`Sminus1' {
			loc stub`r' = subinstr("`stub'","`S'","`r'",1)
			loc index1 = colnumb(`pvector'_s`s',"`stub`r''")
			loc param = `param' + `pvector'_s`s'[1,`index1']
			}
		mat `pvector'_s`s'[1,`index2'] = (0 - `param')
		}
	}
*mat list `pvector'_s`S'

	
** Symmetry

foreach stub in np`S' ynp`S' np`S'z1 {
    forvalues s = 1/`Sminus1' {
	    loc stub2 = subinstr("`stub'","`S'","`s'",1)
		loc index2 = colnumb(`pvector'_s`S',"`stub2'")
	    loc index1 = colnumb(`pvector'_s`s',"`stub'")
		mat `pvector'_s`S'[1,`index2']=`pvector'_s`s'[1,`index1']
		}
	}
*mat list `pvector'_s`S'


** Summing to zero within `S' 

foreach stub in np`S' ynp`S' np`S'z1 {
	loc param = 0
	loc index1 = colnumb(`pvector'_s`S',"`stub'")
    forvalues s = 1/`Sminus1' {
	    loc stub2 = subinstr("`stub'","`S'","`s'",1)
	    loc index2 = colnumb(`pvector'_s`S',"`stub2'")
		loc param = `param' + `pvector'_s`S'[1,`index2']
		}
	mat `pvector'_s`S'[1,`index1'] = (0 - `param')
	}
*mat list `pvector'_s`S'

	
forvalues s = 1/`S' {
	return matrix `pvector'_s`s' = `pvector'_s`s', copy
	}
}
end


********************************************
** Define a program to predict s_hat      **
********************************************

cap program drop s_hat_recover
program define s_hat_recover
args numeq
 
loc neqplus1 = `numeq' + 1
** Predict the first category and recover the last one
tempvar  s_tot
g `s_tot' = 0 /*the dropped category's share will be recovered the other categories' predicted shares*/
forvalues g = 1/`numeq' { 
	predict s`g'_hat, equation(s`g') xb 
	replace `s_tot' = `s_tot' + s`g'_hat
	}
	g s`neqplus1'_hat = 1 - `s_tot'

** Correct out of range predictions
forvalues g = 1/`neqplus1' {
    count if s`g'_hat > 1 /*27 obs for s1*/
	replace s`g'_hat = .999999 if s`g'_hat > 1
	count if s`g'_hat < 0 /*5 obs for s1, 7 obs for s2, 247 obs for s3*/
	replace s`g'_hat = .000001 if s`g'_hat < 0
	
	su s`g'_hat, detail 
	}
end


************************************************
** Define a program to calculate elasticities **
************************************************

cap program drop elasticity_calc
program define elasticity_calc
args numeq 

loc neqplus1 = `numeq' + 1

forvalues g = 1/`neqplus1' {
	g nabla_omega_g`g' = 0
	loc nabla_omega_g`g'_expression 0
	forvalues r = 1/${npowers} {
		loc nabla_omega_g`g'_expression `nabla_omega_g`g'_expression' + params_s`g'[1,colnumb(params_s`g',"y`r'")]*(y1^(`r'-1))*`r'
		replace nabla_omega_g`g' = nabla_omega_g`g' + params_s`g'[1,colnumb(params_s`g',"y`r'")] *(y1^(`r'-1))*`r'
		}
	forvalues k = 1/`neqplus1' {
		replace nabla_omega_g`g' = nabla_omega_g`g' + params_s`g'[1,colnumb(params_s`g',"ynp`k'")]*np`k'
		}
	replace nabla_omega_g`g' = nabla_omega_g`g' + params_s`g'[1,colnumb(params_s`g',"yz1")]*z1 
	}

g pBp_denom = 0
loc pBp 0
forvalues g = 1/`neqplus1' {
	cap drop rowsum
	g rowsum = 0
	loc rowsum 0
	forvalues k = 1/`neqplus1' {
		loc rowsum `rowsum' + params_s`g'[1,colnumb(params_s`g',"ynp`k'")]*np`k'
		replace rowsum = rowsum + params_s`g'[1,colnumb(params_s`g',"ynp`k'")]*np`k'
		}
	loc pBp `pBp' + np`g'*(`rowsum')
*	display `"pBp: `pBp' "'
	replace pBp_denom = pBp_denom + np`g'*(`rowsum')
	}

replace pBp_denom = (1+0.5*pBp)
*su pBp_denom
mat I`neqplus1' = I(`neqplus1')

forvalues g = 1/`neqplus1' { 
    g ds`g'_ex = .
	}

count
loc N = r(N)

quietly {  
	forvalues i = 1/`N' { 
		mat nabla_omega_p = J(`neqplus1',`neqplus1',.)
		mat list nabla_omega_p
		mat part2 = J(`neqplus1',1,.)
		mat list part2
		forvalues g=1/`neqplus1' {
			mat part2[`g',1] = nabla_omega_g`g'[`i']/(pBp_denom[`i'])
			forvalues k = 1/`neqplus1' {
				mat nabla_omega_p[`g',`k'] = nabla_omega_g`g'[`i']*np`k'[`i']
				}
			}
		mat list nabla_omega_p
		mat list part2
		mat ds_ex = inv(I`neqplus1'+(1/pBp_denom[`i'])*nabla_omega_p)*part2
		mat list ds_ex
		forvalues g=1/`neqplus1' {
			replace ds`g'_ex = ds_ex[`g',1] in `i'
			}
		}
	}

cap drop pBP_denom nabla_omega_g* rowsum

forvalues g=1/`neqplus1' { 
	g es`g'_ex = ds`g'_ex * 1/(s`g'_hat)
	g de`g'_ex = exp(y1)*(ds`g'_ex + s`g'_hat) 
	g ee`g'_ex = ds`g'_ex * (1/s`g'_hat) + 1
	}


** Price Elasticities
* Get dsi_eprj epri_eprj
* the hicksian semi- budget share elasticity for good i wrt price j is : dwi_epj = A0_{ij} + B_{ij}*y1 + A1_{ij}*z1
* and the hicksian price elasticity is : eei_epj = dwi_epj/wi - 1_ij + wj 

forvalues g = 1/`neqplus1' {  
	forvalues k = 1/`neqplus1' {
		g ds`g'_epr`k' = params_s`g'[1,colnumb(params_s`g',"np`k'")]+z1*params_s`g'[1,colnumb(params_s`g',"np`k'z1")]+y1*params_s`g'[1,colnumb(params_s`g',"ynp`k'")]
		}
	}
	
forvalues g = 1/`neqplus1' { 
	forvalues k = 1/`neqplus1' {
	    g ee`g'_epr`k' = ds`g'_epr`k'*(1/s`g'_hat) + s`k'_hat
		if `k'==`g' {
			replace ee`g'_epr`k' = ee`g'_epr`k' - 1
			}
		}
	}

forvalues g = 1/`neqplus1' {
	forvalues zs = 1/$nshift {
		g ds`g'_dzs`zs' =  params_s`g'[1,colnumb(params_s`g',"zs`zs'")]
		} 
	}
		
end


****************************************************
** First Load Results and Recover Lost Parameters **
****************************************************


** Main results no instruments

estimates use "${OUT}\reg3_iterated_all_`mainspec'"
mat params = e(b)

loc K = colsof(params)
loc S = $neq + 1 /*last (recovered) equation*/
loc Sminus1 = $neq 

param_recover_easi params $neq

forvalues s = 1/`S' {
    mat params_s`s' = r(params_s`s')
	mat list params_s`s'
	mat params_mean = nullmat(params_mean),params_s`s'
	}


count
loc N = r(N)

loc shat_paramlist "sf_hat"
loc shat_collapselist_sd "sf_hat_se=sf_hat"
loc dsex_paramlist "dsf_ex"
loc dsex_collapselist_sd "dsf_ex_se=dsf_ex"
loc esex_paramlist "esf_ex"
loc esex_collapselist_sd "esf_ex_se=esf_ex"
loc deex_paramlist "def_ex"
loc deex_collapselist_sd "def_ex_se=def_ex"
loc eeex_paramlist "eef_ex"
loc eeex_collapselist_sd "eef_ex_se=eef_ex"
loc dsepr_paramlist ""
loc dsepr_collapselist_sd ""
loc eeepr_paramlist ""
loc eeepr_collapselist_sd ""
loc zs_paramlist ""
loc zs_collapselist_sd ""

forvalues s=1/`S' {   
	loc shat_paramlist `shat_paramlist' s`s'_hat
	loc shat_collapselist_sd `shat_collapselist_sd' s`s'_hat_se=s`s'_hat
	loc dsex_paramlist `dsex_paramlist' ds`s'_ex
	loc dsex_collapselist_sd `dsex_collapselist_sd' ds`s'_ex_se=ds`s'_ex
	loc esex_paramlist `esex_paramlist' es`s'_ex
	loc esex_collapselist_sd `esex_collapselist_sd' es`s'_ex_se=es`s'_ex
	loc deex_paramlist `deex_paramlist' de`s'_ex
	loc deex_collapselist_sd `deex_collapselist_sd' de`s'_ex_se=de`s'_ex
	loc eeex_paramlist `eeex_paramlist' ee`s'_ex
	loc eeex_collapselist_sd `eeex_collapselist_sd' ee`s'_ex_se=ee`s'_ex
	display "did first part share `s'"
	
	forvalues k=1/`S' {
		loc dsepr_paramlist `dsepr_paramlist' ds`s'_epr`k'
		loc dsepr_collapselist_sd `dsepr_collapselist_sd' ds`s'_epr`k'_se=ds`s'_epr`k'
		loc eeepr_paramlist `eeepr_paramlist' ee`s'_epr`k'
		loc eeepr_collapselist_sd `eeepr_collapselist_sd' ee`s'_epr`k'_se=ee`s'_epr`k'
		display "did second part share `s' k `k'"
		}
		
	forvalues zs=1/$nshift {
		loc zs_paramlist `zs_paramlist' ds`s'_dzs`zs'
		loc zs_collapselist_sd `zs_collapselist_sd' ds`s'_dzs`zs'_se=ds`s'_dzs`zs'
		display "did third part share `s' zs `zs'"
		}
	}

	
loc paramlist `dsex_paramlist' `esex_paramlist' `deex_paramlist' `eeex_paramlist' `dsepr_paramlist' `eeepr_paramlist'



************************************
** Post Estimation - Elasticities **
************************************

** First: Expenditure elasticities using recovered parameters  

use "${OUT}\reg3_iterated_data_`mainspec'.dta", clear
g np`S' = 0 /*prices are normalized so last price is zero*/
g x = y1 

forvalues g = 1/`S' {
    replace x = x + p`g'*s`g' /* this is total expenditures (which is log and demeaned) but not deflated*/
	}

foreach pop in urb rur all {	
	xtile qn_`pop' = y1 if `pop'==1, n(5)
	forvalues q = 1/4 {
		su y1 if qn_`pop'==`q'
		loc q`q'_`pop' = r(max)
		}
	}

	
count
loc N = r(N)
loc N : di %9.0fc `N'

su y1
loc y1_mean = r(mean)

forvalues s = 1/`S' {
    su np`s'
	loc np`s'_mean = r(mean)
	}


tempfile estimation_dataset
save `estimation_dataset', replace


** Predict s_hat using program above

s_hat_recover $neq


** Calculate elasticities using program above  

elasticity_calc $neq

su ds1_ex - ee`S'_epr`S', detail
bysort qn_all: su ds1_ex - ee`S'_epr`S'


* Create a weighted average food exp elast
g sf_hat = 0
forvalues f = 1/$nf {
	replace sf_hat = sf_hat + s`f'_hat 
	}

foreach param in ds es de ee { /*arithmetic mean*/
	g `param'f_ex = 0
		forvalues f=1/$nf {
			replace `param'f_ex = `param'f_ex + `param'`f'_ex*s`f'_hat
			}
		replace `param'f_ex = `param'f_ex*(1/sf_hat)
	}
	
tempfile reg3_elast
save `reg3_elast', replace


******************************
** Simulate standard errors **
******************************

** Steps: 
* For each bootstrap iteration b... 
* Load row b of the bootstrapped parameters. 
* Erpost to estimates. 
* Recover parameters. 
* Calculate elasticities.

if `dobootstrap' == 1 {
estimates use "${OUT}\reg3_iterated_all_`mainspec'"
mat params = e(b)
loc eqnames: coleq params
loc colnames: colnames params

loc K = colsof(params)
loc S = $neq + 1 /*last (recovered) equation*/
loc Sminus1 = $neq 

use  "${OUT}\reg3_iterated_params_boot_all_`mainspec'", clear
mkmat _all, matrix(params_boot_all)
mat coleq params_boot_all = `eqnames'
mat colnames params_boot_all = `colnames'

loc matappend "params_s1"
forvalues s=2/`S' {
	loc matappend "`matappend',params_s`s'"
	}

count
loc B = r(N)

forvalues b = 1/`B' {
    display as error "Bootstrap `b' of `B':"
	
	*time stamp
	local today =  "`c(current_date)' `c(current_time)'"
	display as error "Time stamp: `today'"

   	use `estimation_dataset', clear		
	estimates use "${OUT}\reg3_iterated_all_`mainspec'"
 
	mat params = params_boot_all[`b',1..`K']
	erepost b = params
	mat params = e(b)

	param_recover_easi params $neq

	mat params_b_all = `matappend'

	* build a list of all parameters including recovered
	mat params_B_all = nullmat(params_B_all)\params_b_all
	
	s_hat_recover $neq

	elasticity_calc $neq
	
	keep hh_c round s*hat ds1_ex - ds`S'_dzs$nshift 
	
	g sf_hat = 0
	forvalues f = 1/$nf {
		replace sf_hat = sf_hat + s`f'_hat
		}
	
	foreach param in ds es de ee { 
		g `param'f_ex = 0
		forvalues f=1/$nf {
			replace `param'f_ex = `param'f_ex + `param'`f'_ex*s`f'_hat
			}
		replace `param'f_ex = `param'f_ex*(1/sf_hat)
		}
	
	tempfile reg3_elast_boot_`b'
	save `reg3_elast_boot_`b'', replace
	}

	
** Append all of the simulated SEs
clear
forvalues b = 1/`B' {
	append using `reg3_elast_boot_`b''
	}
	

* drop outliers using IQ range method
foreach param in `paramlist' {
	display "`param'"
	count if `param'==.
	su `param', detail
	loc pct25 = r(p25)
	loc pct75 = r(p75)
	replace `param' = . if `param' > `pct75' + 1.5*(`pct75'-`pct25') & `param' != .
	replace `param' = . if `param' < `pct25' - 1.5*(`pct75'-`pct25') & `param' != .
	count if `param'==.
	}
	
** Collapse to the hh-round level 

collapse (sd) `shat_collapselist_sd' `dsex_collapselist_sd' `esex_collapselist_sd' `deex_collapselist_sd' `eeex_collapselist_sd' ///
`dsepr_collapselist_sd' `eeepr_collapselist_sd' `zs_collapselist_sd' ///
, by(hh_c round)


merge 1:1 hh_c round using `reg3_elast', nogen keep(1 3) keepusing(qn* y1 x $zlist $nplist $zslist urb rur all ds1_ex - ds`S'_dzs$nshift sf_hat dsf_ex esf_ex def_ex eef_ex)  

save "${OUT}\reg3_elasticities_`mainspec'", replace

clear
svmat params_B_all
save "${OUT}\reg3_iterated_params_boot_all_recovered_`mainspec'", replace
}


************************
** Elasticity Figures **
************************

** Elasticity figures by quintile

use "${OUT}\reg3_elasticities_`mainspec'", clear
keep if round==3 /*Depict elasticities for most recent round*/
	
loc nbins 10

tempfile reg3_elasticities_bin
save `reg3_elasticities_bin'

foreach pop in urb rur all {
	foreach param in `paramlist' {
		display "`param'"
		
		use `reg3_elasticities_bin', clear

		keep if `pop'==1
		xtile bin = y1, n(`nbins')

		sort bin `param'
		keep bin `param' `param'_se
		
		collapse (median) `param' `param'_se, by(bin)

		g pop = "`pop'"
		keep bin pop `param' `param'_se

		tempfile `param'_medval_`pop'
		save ``param'_medval_`pop'', replace		
		}
	}
	
	
clear
set obs `nbins'
g bin = _n
expand 3
g count = 1
bysort bin : g popc = sum(count) 
g pop = "all" if popc==1
replace pop = "urb" if popc==2
replace pop = "rur" if popc==3
forvalues i=1/`S' {
	g order`i' = (bin-1)*(`S'+1) + `i'	
	}

drop popc count

foreach pop in urb rur all {
	foreach param in `paramlist' {
		merge 1:1 bin pop using ``param'_medval_`pop'', nogen update
		}
	}
	
foreach param in `paramlist' {
	g `param'_sim_lo = `param' - `param'_se*1.96
	g `param'_sim_hi = `param' + `param'_se*1.96		
	}

save "${OUT}\elasticity_summary_byq_`mainspec'", replace




** Plot exp elasticity point estimates with confidence intervals by quintile and population 

loc dslab "Budget share semi-elasticity"
loc dsyline 0
loc ymin_ds -0.3
loc ymax_ds 0.35
loc yskip_ds 0.1
loc eslab "Budget share elasticity"
loc esyline 0
loc ymin_es -1.5
loc ymax_es 2.5
loc yskip_es 0.5
loc delab "Expenditure semi-elasticity"
loc deyline 0
loc ymin_de -0.5
loc ymax_de 5
loc yskip_de 1
loc eelab "Expenditure elasticity"
loc eeyline 1
loc ymin_ee -.5
loc ymax_ee 3.5
loc yskip_ee 0.5

loc xlab ""
loc o = (`S'+1)/2
forvalues b=1/`nbins' {
    loc xlab `xlab' `o' "bin `b'"
	loc o = `o'+(`S'+1)
	}
display `"The label is: `xlab'"'

foreach param in ds es de ee {
	foreach pop in urb rur all {
		graph twoway 	(rcap `param'1_ex_sim_hi `param'1_ex_sim_lo order1, color(red%20)) ///
						(scatter `param'1_ex order1, color(red%20))  ///
						(rcap `param'2_ex_sim_hi `param'2_ex_sim_lo order2, color(green%50)) ///
						(scatter `param'2_ex order2, color(green%50)) ///
						(rcap `param'3_ex_sim_hi `param'3_ex_sim_lo order3, color(blue%75)) ///
						(scatter `param'3_ex order3, color(blue%75)) ///
						if pop=="`pop'" ///
						, legend(order(2 "Food" 4 "Goods" 6 "Services") rows(1)) title("``param'lab' w.r.t. to total exp., `pop' pop") xlabel(`xlab', noticks) xtitle("Log real expenditures") yline(``param'yline') ylabel(`ymin_`param'' (`yskip_`param'') `ymax_`param'', nogrid) yscale(range(`ymin_`param'' `ymax_`param''))
		graph export "${FIG}\Elast_byq_`param'ee_`pop'_`mainspec'.pdf", replace  
		}
	}
	

** Plot price elasticity point estimates with confidence intervals by quintile and population 

loc dslab "Budget share semi-elasticity"
loc dsyline 0
loc ymin_ds_epr -0.5
loc ymax_ds_epr 0.5
loc yskip_ds_epr 0.25
loc eelab "Hicksian price elasticity"
loc eeyline 0
loc ymin_ee_epr -3
loc ymax_ee_epr 3
loc yskip_ee_epr 0.5

	
foreach param in ds ee {
	foreach pop in urb rur all {
		forvalues g = 1/`S' {
			graph twoway 	(rcap `param'1_epr`g'_sim_hi `param'1_epr`g'_sim_lo order1, color(red%20)) ///
							(scatter `param'1_epr`g' order1, color(red%20))  ///
							(rcap `param'2_epr`g'_sim_hi `param'2_epr`g'_sim_lo order2, color(green%50)) ///
							(scatter `param'2_epr`g' order2, color(green%50)) ///
							(rcap `param'3_epr`g'_sim_hi `param'3_epr`g'_sim_lo order3, color(blue%90)) ///
							(scatter `param'3_epr`g' order3, color(blue%90)) ///
							if pop=="`pop'" ///
							, legend(order(2 "Food" 4 "Goods" 6 "Services") rows(1)) title("``param'lab' with respect to `g`g'' price, `pop' pop") xlabel(`xlab', noticks) xtitle("Log real expenditures") yline(``param'yline') ylabel(`ymin_`param'_epr' (`yskip_`param'_epr') `ymax_`param'_epr', nogrid) yscale(range(`ymin_`param'_epr' `ymax_`param'_epr'))
			graph export "${FIG}\PriceElast_byq_`param'_epr_`g`g''_`pop'_`mainspec'.pdf", replace  
			}
		}
	}
	

** Elasticity lpoly regresssions

loc dsexlab "Budget share semi-elasticity"
use "${OUT}\reg3_elasticities_`mainspec'",  clear
keep if round==3 

su y1, detail
loc y1_min = r(p1)
display "`y1_min'"
loc y1_max = r(p99)
display "`y1_max'"

foreach pop in urb rur all {	
	forvalues q = 1/4 {
		su y1 if qn_`pop'==`q'
		loc q`q'_`pop' = r(max)
		}		
	}
	
kdensity y1 if y1 >= `y1_min' & y1 <= `y1_max', nograph generate(grid_y1 predd_y1)

tempfile reg3_lpoly
save `reg3_lpoly', replace

loc plot1 ds1_ex ds2_ex ds3_ex 
loc plotlab1 "reg3_lpoly_dsg_ex"
loc plottitle1 "Budget share semi-elasticity"
loc ymin_plot1 = -0.5
loc ymax_plot1 = 0.6
loc yline_plot1 = 0
loc yskip_plot1 = 0.1

loc plot2 ee1_ex ee2_ex ee3_ex 
loc plotlab2 "reg3_lpoly_eeg_ex"
loc plottitle2 "Total expenditure elasticity"
loc ymin_plot2 = -4
loc ymax_plot2 = 4
loc yline_plot2 = 1
loc yskip_plot2 = 1

loc plot3 ee1_epr1 ee2_epr1 ee3_epr1
loc plotlab3 "reg3_lpoly_eeg_epr1"
loc plottitle3 "Elasticity of demand wrt ``g1'lab' price"
loc ymin_plot3 = -1.5
loc ymax_plot3 = 1.5
loc yline_plot3 = 0
loc yskip_plot3 = 0.5

loc plot4 ee1_epr2 ee2_epr2 ee3_epr2
loc plotlab4 "reg3_lpoly_eeg_epr2"
loc plottitle4 "Elasticity of demand wrt ``g2'lab' price"
loc ymin_plot4 = -1.5
loc ymax_plot4 = 1.5
loc yline_plot4 = 0
loc yskip_plot4 = 0.5

loc plot5 ee1_epr3 ee2_epr3 ee3_epr3
loc plotlab5 "reg3_lpoly_eeg_epr3"
loc plottitle5 "Elasticity of demand wrt ``g3'lab' price"
loc ymin_plot5 = -1.5
loc ymax_plot5 = 1.5
loc yline_plot5 = 0
loc yskip_plot5 = 0.5

loc c1 red%15
loc c2 green%50 
loc c3 blue%80 

forvalues plot = 1/1 {
	loc paramlist `plot`plot''
	
	cap drop lpoly_* 
	
	foreach pop in urb rur all {	
		use `reg3_lpoly', clear
		
		loc command ""
		
		loc i = 1

		foreach param in `paramlist' {
			cap drop `param'_lo `param'_hi
			g `param'_lo = `param' - `param'_se*1.96 if `pop'==1 & y1 >= `y1_min' & y1 <= `y1_max'
			g `param'_hi = `param' + `param'_se*1.96 if `pop'==1 & y1 >= `y1_min' & y1 <= `y1_max'

			lpoly `param' y1 if 	`pop'==1 & 	y1 >= `y1_min' & 	y1 <= `y1_max', nograph degree(0) at(grid_y1) generate(lpoly_`plot'_`param'_`pop')
			lpoly `param'_lo y1 if 	`pop'==1 & 	y1 >= `y1_min' & 	y1 <= `y1_max', nograph degree(0) at(grid_y1) generate(lpoly_`plot'_`param'_`pop'_lo)
			lpoly `param'_hi y1 if 	`pop'==1 & 	y1 >= `y1_min' & 	y1 <= `y1_max', nograph degree(0) at(grid_y1) generate(lpoly_`plot'_`param'_`pop'_hi)
			
			loc command `command' (line lpoly_`plot'_`param'_`pop' grid_y1, lpattern(solid) lcolor(`c`i'') lwidth(medthick)) (line lpoly_`plot'_`param'_`pop'_lo grid_y1, lpattern(dash) lcolor(`c`i'') lwidth(vthin)) (line lpoly_`plot'_`param'_`pop'_hi grid_y1, lpattern(dash) lcolor(`c`i'') lwidth(vthin))
			
			loc i = `i'+1
			}
		
		graph twoway `command' ///
			||, ///
			xscale(range(`y1_min' `y1_max')) xline(`q1_`pop'', lwidth(vthin) lpattern(dash) lcolor(gray)) xline(`q2_`pop'', lwidth(vthin) lpattern(dash) lcolor(gray)) xline(`q3_`pop'', lwidth(vthin) lpattern(dash) lcolor(gray)) xline(`q4_`pop'', lwidth(vthin) lpattern(dash) lcolor(gray))  ///
			yscale(range(`ymin_plot`plot'' `ymax_plot`plot'')) ylabel(`ymin_plot`plot'' (`yskip_plot`plot'') `ymax_plot`plot'', nogrid) yline(`yline_plot`plot'', lcolor(gray) lwidth(medthick)) ///
			title("`plottitle`plot'' (`pop')") ///
			ytitle("") ///
			legend(rows(1) order(1 "``g1'lab'" 4 "``g2'lab'" 7 "``g3'lab'" ))
		
		graph export "${FIG}\reg3_lpoly_`plotlab`plot''_`pop'_`mainspec'.pdf", replace 
		}
	}

	
** Make table of all parameters 

loc colnames : colnames params_mean
loc coleqnames : coleq params_mean
loc paramlist : colnames params_s1

use "${OUT}\reg3_iterated_params_boot_all_recovered_`mainspec'", clear
mkmat _all, matrix(params_B_all)
count
loc B = r(N)

mat colnames params_B_all = `colnames'
mat coleq params_B_all = `eqnames'

loc K = colsof(params_mean)
mat vce_tot = J(`K',`K',0)
forvalues b=1/`B' {
	mat vce_b = (params_B_all[`b',1...]-params_mean)'*(params_B_all[`b',1...]-params_mean)
	mat vce_tot = vce_tot + vce_b
	}

mat vce_boot=J(`K',`K',.)
forvalues i = 1/`K' {
	forvalues j = 1/`K' {
		mat vce_boot[`i',`j']=(1/(`B'-1))*vce_tot[`i',`j']
		}
	}

mat params_se = J(1,`K',.)
mat colnames params_se = `colnames'
mat coleq params_se = `coleqnames'

forvalues i = 1/`K' {
	mat params_se[1,`i'] = sqrt(vce_boot[`i',`i'])
	}
	

* specify variable labels:
loc y1_lab "\( Y \)"
forvalues j = 2/$npowers {
	loc y`j'_lab "\( Y^`j' \) "
	}
loc z1_lab "Urban"
loc yz1_lab "Y x urban"
forvalues g=1/`S' {
	loc np`g'_lab "``g`g''lab' price"
	loc ynp`g'_lab "Y x ``g`g''lab' price"
	loc np`g'z1_lab "``g`g''lab' price x urban"
	}
	
loc zs1_lab "HH size (adult equiv)"
loc zs2_lab "Consumptive asset index" 
loc zs3_lab "Educ of hh head (yrs)"
loc zs4_lab "Max educ by hh member (yrs)"
loc zs5_lab "HH dependency ratio (share)"
loc zs6_lab "HH part. farming"
loc zs7_lab "HH part. non-farm self empl"
loc zs8_lab	"HH part. wage labor"
loc zs18_lab "Food price inst 1" 
loc zs19_lab "Food price inst 2"
loc zs20_lab "Food price inst 3"
loc zs21_lab "Goods price inst 1"
loc zs22_lab "Goods price inst 2"
loc zs23_lab "Goods price inst 3"
loc zs24_lab "Services price inst 1"
loc zs25_lab "Services price inst 2"
loc zs26_lab "Services price inst 3"

display "`colnames'"

cap texdoc close 
texdoc init "${FIG}\Parameters_easi_reg3_`mainspec'", replace force

display "N: `N'"
loc cv05 = 1.960 /*invttail(`N',0.05/2)*/
loc cv01 = 2.576 /*invttail(`N',0.01/2)*/
loc cv001 = 3.292 /*invttail(`N',0.001/2)*/

foreach var in `paramlist' {
    display "Var: `var'"
	
	if !(("`var'"=="np`S'") | ("`var'"=="ynp`S'") | ("`var'"=="np`S'z1") | ("`var'"=="zs0") | ("`var'"=="zs9") | ("`var'"=="zs10") | ("`var'"=="zs11") | ("`var'"=="zs12") | ("`var'"=="zs13") | ("`var'"=="zs14") | ("`var'"=="zs15") | ("`var'"=="zs16") | ("`var'"=="zs17") | ("`var'"=="crep1") | ("`var'"=="crep2") | ("`var'"=="crep3") | ("`var'"=="crep4") | ("`var'"=="crep5") | ("`var'"=="crep6") | ("`var'"=="crepy1") | ("`var'"=="crepy2") | ("`var'"=="crepy3") | ("`var'"=="crepy4") | ("`var'"=="crepy5") | ("`var'"=="crepy6")) {

		display "Var: `var'; Including in table."		
		
		forvalues s = 1/`S' {
			loc mean_var_s`s' : di %9.3f params_mean[1,"s`s':`var'"]
			loc se_var_s`s' : di %9.3f params_se[1,"s`s':`var'"]
			loc sig_var_s`s' = (abs(`mean_var_s`s''/`se_var_s`s'')>`cv05') + (abs(`mean_var_s`s''/`se_var_s`s'')>`cv01') + (abs(`mean_var_s`s''/`se_var_s`s'')>`cv001')
			loc stars_var_s`s' ""
			if `sig_var_s`s''==1 {
				loc stars_var_s`s' "*"
				}
				else if `sig_var_s`s''==2 {
					loc stars_var_s`s' "**"
					}
					else if `sig_var_s`s''==3 {
						loc stars_var_s`s' "***"
						}
			}
			
		tex ``var'_lab' 	& `mean_var_s1'`stars_var_s1' & (`se_var_s1') & `mean_var_s2'`stars_var_s2' & (`se_var_s2') & `mean_var_s3'`stars_var_s3' & (`se_var_s3')  \\
		}
	}

loc var zs0
forvalues s = 1/`S' {
	loc mean_var_s`s' : di %9.3f params_mean[1,"s`s':`var'"]
	loc se_var_s`s' : di %9.3f params_se[1,"s`s':`var'"]
	loc sig_var_s`s' = (abs(`mean_var_s`s''/`se_var_s`s'')>`cv05') + (abs(`mean_var_s`s''/`se_var_s`s'')>`cv01') + (abs(`mean_var_s`s''/`se_var_s`s'')>`cv001')
	loc stars_var_s`s' ""
	if `sig_var_s`s''==1 {
		loc stars_var_s`s' "*"
		}
		else if `sig_var_s`s''==2 {
			loc stars_var_s`s' "**"
			}
			else if `sig_var_s`s''==3 {
				loc stars_var_s`s' "***"
				}
	}	
tex Intercept 	& `mean_var_s1'`stars_var_s1' & (`se_var_s1') & `mean_var_s2'`stars_var_s1' & (`se_var_s2') & `mean_var_s3'`stars_var_s1' & (`se_var_s3') \\

tex \hline
tex N: & \multicolumn{6}{c}{`N'} \\ 
tex Fixed effects: & \multicolumn{6}{c}{Macro-region, survey round} \\ 
tex Correlated random effects: & \multicolumn{6}{c}{Household level} \\
tex & \multicolumn{6}{c}{(Prices and price-total expenditure interactions)} \\ \hline 

texdoc close


** Make table of budget share semi-elasticities

use "${OUT}\reg3_elasticities_`mainspec'", clear

texdoc init "${FIG}\AMEs_easi_reg3_`mainspec'", replace force

foreach pop in all rur urb {
	forvalues s = 1/`S' {
		su ds`s'_ex if `pop'==1
		loc y1_ame_s`s'_`pop'_mean : di %9.3f r(mean)
		su ds`s'_ex_se if `pop'==1
		loc y1_ame_s`s'_`pop'_se : di %9.3f r(mean)

		loc sig_y1_ame_s`s'_`pop' = (abs(`y1_ame_s`s'_`pop'_mean'/`y1_ame_s`s'_`pop'_se')>`cv05') + (abs(`y1_ame_s`s'_`pop'_mean'/`y1_ame_s`s'_`pop'_se')>`cv01') + (abs(`y1_ame_s`s'_`pop'_mean'/`y1_ame_s`s'_`pop'_se')>`cv001')
			loc stars_y1_ame_s`s'_`pop' ""
			if `sig_y1_ame_s`s'_`pop''==1 {
				loc stars_y1_ame_s`s'_`pop' "*"
				}
				else if `sig_y1_ame_s`s'_`pop''==2 {
					loc stars_y1_ame_s`s'_`pop' "**"
					}
					else if `sig_y1_ame_s`s'_`pop''==3 {
						loc stars_y1_ame_s`s'_`pop' "***"
						}

		}
	
	forvalues t = 1/`S' {
		forvalues s = 1/`S' {	
			su ds`s'_epr`t' if `pop'==1
			loc np`t'_ame_s`s'_`pop'_mean : di %9.3f r(mean)
			su ds`s'_epr`t'_se if `pop'==1
			loc np`t'_ame_s`s'_`pop'_se : di %9.3f r(mean)
			
			loc sig_np`t'_ame_s`s'_`pop' = (abs(`np`t'_ame_s`s'_`pop'_mean'/`np`t'_ame_s`s'_`pop'_se')>`cv05') + (abs(`np`t'_ame_s`s'_`pop'_mean'/`np`t'_ame_s`s'_`pop'_se')>`cv01') + (abs(`np`t'_ame_s`s'_`pop'_mean'/`np`t'_ame_s`s'_`pop'_se')>`cv001')
			loc stars_np`t'_ame_s`s'_`pop' ""
			if `sig_np`t'_ame_s`s'_`pop''==1 {
				loc stars_np`t'_ame_s`s'_`pop' "*"
				}
				else if `sig_np`t'_ame_s`s'_`pop''==2 {
					loc stars_np`t'_ame_s`s'_`pop' "**"
					}
					else if `sig_np`t'_ame_s`s'_`pop''==3 {
						loc stars_np`t'_ame_s`s'_`pop' "***"
						}
			}
		}
		
	foreach z in 1 2 3 4 5 6 7 8  {
	    forvalues s = 1/`S' {
		    su ds`s'_dzs`z' if `pop'==1 
			loc zs`z'_ame_s`s'_`pop'_mean : di %9.3f r(mean)
			su ds`s'_dzs`z'_se if `pop'==1
			loc zs`z'_ame_s`s'_`pop'_se : di %9.3f r(mean)
			
			loc sig_zs`z'_ame_s`s'_`pop' = (abs(`zs`z'_ame_s`s'_`pop'_mean'/`zs`z'_ame_s`s'_`pop'_se')>`cv05') + (abs(`zs`z'_ame_s`s'_`pop'_mean'/`zs`z'_ame_s`s'_`pop'_se')>`cv01') + (abs(`zs`z'_ame_s`s'_`pop'_mean'/`zs`z'_ame_s`s'_`pop'_se')>`cv001')
			loc stars_zs`z'_ame_s`s'_`pop' ""
			if `sig_zs`z'_ame_s`s'_`pop''==1 {
				loc stars_zs`z'_ame_s`s'_`pop' "*"
				}
				else if `sig_zs`z'_ame_s`s'_`pop''==2 {
					loc stars_zs`z'_ame_s`s'_`pop' "**"
					}
					else if `sig_zs`z'_ame_s`s'_`pop''==3 {
						loc stars_zs`z'_ame_s`s'_`pop' "***"
						}
			}
		}
	loc nshiftcoreplus1 = ${nshift_core} + 1
	loc nshiftplusexogprices = ${nshift}
	display "nshiftcoreplus1: `nshiftcoreplus1' nshiftplusexogprices: `nshiftplusexogprices'"
	cap su ds1_dzs`nshiftcoreplus1'
	if _rc==0 {
		forvalues z = `nshiftcoreplus1'/`nshiftplusexogprices' { 
			forvalues s = 1/`S' {
				su ds`s'_dzs`z' if `pop'==1 
				loc zs`z'_ame_s`s'_`pop'_mean : di %9.3f r(mean)
				su ds`s'_dzs`z'_se if `pop'==1
				loc zs`z'_ame_s`s'_`pop'_se : di %9.3f r(mean)
				
				loc sig_zs`z'_ame_s`s'_`pop' = (abs(`zs`z'_ame_s`s'_`pop'_mean'/`zs`z'_ame_s`s'_`pop'_se')>`cv05') + (abs(`zs`z'_ame_s`s'_`pop'_mean'/`zs`z'_ame_s`s'_`pop'_se')>`cv01') + (abs(`zs`z'_ame_s`s'_`pop'_mean'/`zs`z'_ame_s`s'_`pop'_se')>`cv001')
			loc stars_zs`z'_ame_s`s'_`pop' ""
			if `sig_zs`z'_ame_s`s'_`pop''==1 {
				loc stars_zs`z'_ame_s`s'_`pop' "*"
				}
				else if `sig_zs`z'_ame_s`s'_`pop''==2 {
					loc stars_zs`z'_ame_s`s'_`pop' "**"
					}
					else if `sig_zs`z'_ame_s`s'_`pop''==3 {
						loc stars_zs`z'_ame_s`s'_`pop' "***"
						}
				} 
			}
		}
		
	count if `pop'==1
	loc N_`pop' : di %9.0fc r(N)
	}
	
loc pop_label_all "Full Sample"
loc pop_label_rur "Rural Observations"
loc pop_label_urb "Urban Observations"
loc ncol_3g 8
loc ncol_6g 16
loc ncolminus2 = `ncol_3g' - 2

foreach pop in rur urb {
	tex \multicolumn{`ncol_3g'}{l}{`pop_label_`pop'' - Expenditure and Price average marginal effects} \\
	tex & Log real expenditures  & `y1_ame_s1_`pop'_mean'`stars_y1_ame_s1_`pop'' & (`y1_ame_s1_`pop'_se')  & `y1_ame_s2_`pop'_mean'`stars_y1_ame_s2_`pop''  & (`y1_ame_s2_`pop'_se')  & `y1_ame_s3_`pop'_mean'`stars_y1_ame_s3_`pop'' & (`y1_ame_s3_`pop'_se')   \\

	forvalues t = 1/`S' {
		tex  & ``g`t''lab' price & `np`t'_ame_s1_`pop'_mean'`stars_np`t'_ame_s1_`pop'' & (`np`t'_ame_s1_`pop'_se') & `np`t'_ame_s2_`pop'_mean'`stars_np`t'_ame_s2_`pop'' & (`np`t'_ame_s2_`pop'_se') & `np`t'_ame_s3_`pop'_mean'`stars_np`t'_ame_s3_`pop'' & (`np`t'_ame_s3_`pop'_se')  \\
		}
	tex \cline{3-`ncol_3g'}
	}


loc pop all
tex \multicolumn{`ncol_3g'}{l}{`pop_label_`pop'' - Demand Shifter average marginal effects} \\
foreach z in 1 2 3 4 5 6 7 8 /*18 19*/  {
	tex & `zs`z'_lab' & `zs`z'_ame_s1_`pop'_mean'`stars_zs`z'_ame_s1_`pop''  & (`zs`z'_ame_s1_`pop'_se') & `zs`z'_ame_s2_`pop'_mean'`stars_zs`z'_ame_s2_`pop'' & (`zs`z'_ame_s2_`pop'_se') & `zs`z'_ame_s3_`pop'_mean'`stars_zs`z'_ame_s3_`pop'' & (`zs`z'_ame_s3_`pop'_se')  \\
	}

cap su ds1_dzs`nshiftcoreplus1'
if _rc==0 {
	forvalues z = `nshiftcoreplus1'/`nshiftplusexogprices' { 
		tex & `zs`z'_lab' & `zs`z'_ame_s1_`pop'_mean'`stars_zs`z'_ame_s1_`pop''  & (`zs`z'_ame_s1_`pop'_se') & `zs`z'_ame_s2_`pop'_mean'`stars_zs`z'_ame_s2_`pop'' & (`zs`z'_ame_s2_`pop'_se') & `zs`z'_ame_s3_`pop'_mean'`stars_zs`z'_ame_s3_`pop'' & (`zs`z'_ame_s3_`pop'_se')  \\
		}
	}
tex \hline
	
	
tex \multicolumn{2}{l}{N (`pop_label_rur'):} & \multicolumn{`ncolminus2'}{c}{`N_rur'} \\
tex \multicolumn{2}{l}{N (`pop_label_urb'):} & \multicolumn{`ncolminus2'}{c}{`N_urb'} \\
tex \multicolumn{2}{l}{Fixed effects:} & \multicolumn{`ncolminus2'}{c}{Macro-region, survey round} \\ 
tex \multicolumn{2}{l}{Correlated random effects:} & \multicolumn{`ncolminus2'}{c}{Household level} \\ 
tex & & \multicolumn{`ncolminus2'}{c}{(Prices and price-total expenditure interactions)} \\ \hline 

texdoc close


log close